1 Loading Packages + Data

## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//Rtmp2PJAwp/downloaded_packages
load(here("jk_code", "JK_cleanMD.rds"))

2 Removing mt, Rp, Gm genes

SO4 <- subset(SO4, features = grep("^mt-|^Rp|^Gm", rownames(SO4), invert = TRUE, value = TRUE))
ElbowPlot(SO4)

SO4 <- SCTransform(SO4) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:27) %>%
    FindClusters(resolution = 0.16) %>%
    RunUMAP(dims = 1:27)
## Running SCTransform on assay: RNA
## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 15822 by 11431
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 286 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 15822 genes
## Computing corrected count matrix for 15822 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 15.90071 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Set default assay to SCT
## PC_ 1 
## Positive:  Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Fabp3, Klk1, Prdx5, Foxq1, Cox6c 
##     Wfdc15b, Ly6a, Wfdc2, Krt7, Slc25a5, Ckb, Cox5b, Atp5g1, Cldn19, Cox4i1 
##     Ggt1, Atp1a1, Ndufa4, Cox8a, Atp5b, Chchd10, Gadd45g, Cox6b1, Atp1b1, Cox7b 
## Negative:  Pappa2, Zfand5, Aard, Robo2, Pde10a, Wwc2, Itga4, Nadk2, Nos1, S100g 
##     Ramp3, Neat1, Sgms2, Ptgs2, Col4a4, Irx1, Col4a3, Tmem158, Mir6236, Ranbp3l 
##     Itprid2, Bmp3, Cdkn1c, Camk2d, Sdc4, Mcub, Dctd, Etnk1, Srrm2, Peg3 
## PC_ 2 
## Positive:  Fos, Junb, Jun, Hspa1b, Hspa1a, Btg2, Egr1, Zfp36, Atf3, Ier2 
##     Fosb, Socs3, Klf6, Klf2, Jund, Dnajb1, Dusp1, Tsc22d1, Gadd45g, Rhob 
##     H3f3b, Gadd45b, Cebpd, Ubc, Mt1, Mt2, Actb, H2bc4, H1f2, Klf4 
## Negative:  Mir6236, Pappa2, CT010467.1, Egf, Slc12a1, Umod, Etnk1, Wnk1, Sfrp1, Atp1b1 
##     Nme7, Robo2, Sptbn1, Col4a4, Hsp90b1, Sec14l1, Pde10a, Zbtb20, Itprid2, Dst 
##     Mal, App, Lars2, Wwc2, Col4a3, Kcnq1ot1, Atrx, Utrn, Tfcp2l1, Pou3f3 
## PC_ 3 
## Positive:  Fth1, Ubb, Ftl1, Ldhb, Car15, Fxyd2, Prdx1, Cd63, Cd9, Eif1 
##     Mpc2, Mgst1, Mt1, Clu, Bsg, Tmem213, Aldoa, Mdh1, Ramp3, Itm2b 
##     Spp1, S100a1, Selenow, Tmem59, Tmem176b, Wfdc2, Tmbim6, Tspo, Atpif1, Ppp1r1a 
## Negative:  Mir6236, Egf, CT010467.1, Umod, Neat1, Tmem52b, Nme7, Fos, Malat1, Slc12a1 
##     Kcnq1ot1, Jun, Junb, Lars2, Etnk1, Dst, Egr1, Wnk1, Slc5a3, Sult1d1 
##     Foxq1, Atrx, Atp1b1, Pnisr, Fosb, Atf3, Ivns1abp, Zfp36, Btg2, Syne2 
## PC_ 4 
## Positive:  Pappa2, Aard, Tmem52b, Umod, Egf, Sult1d1, Tmem158, Foxq1, Mcub, Ptgs2 
##     Dctd, Wwc2, Iyd, Car15, Ramp3, Ptprz1, Hsp90b1, Cd9, S100g, Wnk1 
##     Cdkn1c, Defb42, 1700028P14Rik, Pth1r, 5330417C22Rik, Tmsb4x, Tagln2, Ctsc, Clu, Bmp2 
## Negative:  Mt1, Apoe, Mt2, Aebp1, Sostdc1, Gpx6, Fxyd2, Ptger3, Neat1, Ckb 
##     Fgf9, Egfl6, Ivns1abp, Mfsd4a, Defb1, Car4, Plet1, CT010467.1, Fkbp11, Chchd10 
##     Atp5md, 2900052N01Rik, Cox6c, Igfbp5, Atp5k, Tmem213, Chka, Mgst3, Avpr1a, Atp1a1 
## PC_ 5 
## Positive:  Hspa1a, Hspa1b, CT010467.1, Pappa2, Klk1, Car15, Cldn10, Fth1, Hspa8, Jun 
##     Lars2, Fau, Mir6236, Itm2b, Wfdc2, Uba52-ps, Aard, Ftl1, Ptger3, Eef1a1 
##     Pik3r1, Egf, Sfrp1, Mal, Gpx4, Tspo, Atp1a1, Tmem176a, Id3, Hsp90aa1 
## Negative:  S100g, Actb, Tmem52b, Sdc4, Tmsb10, Ppia, Abhd2, Uroc1, Egr1, Serf2 
##     Slc39a1, Igfbp7, Atf3, Ndufb1-ps, Atp5md, Alkbh5, Cebpd, Cox6c, Ramp3, Ndufa3 
##     Gnb1, Atp5e, Foxq1, Atp5k, Ldhb-ps, Fam107a, Kdm6b, Rbm47, Atp5mpl, Wfdc15b
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11431
## Number of edges: 375699
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8945
## Number of communities: 5
## Elapsed time: 0 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 16:05:11 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:05:11 Read 11431 rows and found 27 numeric columns
## 16:05:11 Using Annoy for neighbor search, n_neighbors = 30
## 16:05:11 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:05:12 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//Rtmp2PJAwp/fileee8534648a7c
## 16:05:12 Searching Annoy index using 1 thread, search_k = 3000
## 16:05:13 Annoy recall = 100%
## 16:05:14 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 16:05:14 Initializing from normalized Laplacian + noise (using RSpectra)
## 16:05:14 Commencing optimization for 200 epochs, with 503018 positive edges
## 16:05:14 Using rng type: pcg
## 16:05:17 Optimization finished
DimPlot(SO4)

DimPlot(SO4,split.by = "sample")

DimPlot(SO4,split.by = "treatment")

SO.markers <- FindAllMarkers(SO4, only.pos = TRUE)
## Calculating cluster 0
## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10

DoHeatmap(SO4, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Phkg1,
## Rtp4, Oasl2, Rsad2, Oasl1, Pvalb, Zfp423, Speer4e, Insyn2a

DimPlot(SO4)

DimPlot(SO4,split.by = "sample")

SO4@meta.data <- SO4@meta.data %>% 
  mutate(subclass_MD = dplyr::case_when(
    seurat_clusters == 0  ~ "type_1",
    seurat_clusters == 1  ~ "type_2",
    seurat_clusters == 2  ~ "type_3",
    seurat_clusters == 3  ~ "type_4",
    seurat_clusters == 4  ~ "type_5",
  

  ))

SO4@meta.data$subclass_MD <- factor(SO4@meta.data$subclass_MD , levels = c("type_1", "type_2", "type_3", "type_4","type_5"))

Idents(SO4) <- SO4@meta.data$subclass_MD

DimPlot(object = SO4, reduction = "umap", group.by = "subclass_MD", label = TRUE)

DimPlot(object = SO4, reduction = "umap", label = TRUE)

Idents(SO4) <- "subclass_MD"

DimPlot(SO4,split.by = "sample",group.by = "seurat_clusters")

DoHeatmap(SO4, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Phkg1,
## Rtp4, Oasl2, Rsad2, Oasl1, Pvalb, Zfp423, Speer4e, Insyn2a

DimPlot(SO4)

markerstoplot <- c("Pappa2", "Egf", "Jun","S100g","Cxcl10")

DotPlot(SO4,
features = markerstoplot,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()

3 Visualizing the DEG’s in each cluster

VlnPlot(SO4,"S100g")
## Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

VlnPlot(SO4,"Cxcl10")

VlnPlot(SO4,"Fos")

VlnPlot(SO4,"Egf")

VlnPlot(SO4,"Pappa2")

save(SO4, file = here("jk_code", "JK_remove_mtrpgm.rds"))